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Abstract 

In order to come closer to a realistic model of high-energy collisions, we simulate SU(2) lattice 
gauge theory under fluctuating temperature. The fluctuations are Euler-Gamma distributed, lead- 
ing to a canonical state maximizing the Renyi and Tsallis entropy formulas. This choice conforms 
to the multiplicity distributions leading to the KNO scaling in high energy experimental spectra. 
We test the random lattice spacing method numerically by investigating the Polyakov Loop ex- 
pectation value, known to be a good order parameter for the confinement - deconfinement phase 
transition in ordinary canonical Monte Carlo methods. The critical coupling (and presumably 
the temperature) move with the width parameter of the inverse temperature fluctuations towards 
higher values. 
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I. INTRODUCTION 



Lattice gauge theory is up to now the only successful nonpertubative numerical approach 
to solve physical problems related to the strong interaction. Among the most reknown recent 
results the prediction of a critical endpoint of the phase transition in QCD became in the 
forefront of research {l -4 1. Also a large scale experimental program, FAIR at GSI, has been 
initiated, among other goals for studying the interface between quark- and hadronic matter 
in the CBM experiment j^]. Accelerator experiments, however, do not have a control on 
thermodynamically relevant parameters, like the temperature and pressure, to such a degree 
that these could be regarded as having a sharp and constant value during the evolution of 
the strongly interacting matter. Lattice theoretical simulations on the other hand assume a 
fixed value for the temperature. 

Our aim with the study presented in this paper is to move towards a more flexible 
scheme: we treat temperature as a random variable, defined not only by its expectation value, 
but also by a width. In fact the thermodynamically consequent approach to this problem 
requires that the inverse temperature, — 1/ ksT, occurring also as a Lagrange multiplier 



for the fixed energy constraint by maximizing t 
randomized. Such a superstatistical method [6 



re entropy, is fixed on the average and then 



11] is in accord with recent findings on non- 



extensive thermodynamics, where the canonical energy distribution is not-exponential, but 
rather shows an experimentally observed power-law tail 12Hl7|. 

In this paper we review basic thermodynamic arguments to relate the temperature to the 
parameters of a statistical power-law tailed, canonical energy distribution. Following this the 
superstatistical method is presented, in particular its realization strategy for lattice Monte 
Carlo simulations. We choose to randomize the timelike to spacelike lattice spacing ratio, 
6 = a t /a s . The most important first task is to check the deconfinement phase transition by 
observing the Polyakov loop expectation value. These results are presented and discussed. 
As a main consequence we predict that the deconfinement transition temperature is likely 
to be higher than determined by fixed-T lattice calculations so far. 
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II. THERMODYNAMIC AL BACKGROUND 



Based on arguments regarding the compatibility of general 
entropy and energy of composed thermodynamical systems 



composition rules for the total 



3 



in an extended canonical 



thermal equilibrium problem the absolute temperature is given by 

= l/T = dL(S)/dL(E), (1) 

with L(S) and L(E) being the additive formal logarithms of the respective composition 
formulas. The formal logarithm maps a general composition law, say Si 2 = Si © S2, to the 
addition by L{S\2) = L(Si) + L(S2). This construction leads us to maximize L(S) — (3L(E) 
when looking for canonical energy distributions [16]. The probability distribution, Wi, of 
states with energy Ei in equilibrium maximizes the formal logarithm of the non-extensive 
entropy formula with constraints on the average value of the also non-additive energy and 
the probability normalization: 

L(S) [wi] - f3^w i L(E i ) - a^Wi = max. (2) 

i i 

Here (3 and a are Lagrange multipliers and it can be proven that /3 = 1/T is related to 
the thermodynamically valid temperature according to the zeroth law of thermodynamics. 
Choosing the next to simplest composition formula to the addition, supplemented with a 
leading second order correction, 

S12 = Si + S 2 + aSiS 2 , (3) 

the additive formal logarithm function is given by 

L(S) = -ln(l + aS). (4) 
a 



This way I/(Si 2 ) = £(Si) + L{S 2 ), indeed. By using the Tsallis entropy formula 181422] 



a 

this formal logarithm turns out to be the Renyi entropy [23 



s = iE(V~ a -^), (5) 



24] 



L(S) = iln5>i-<\ (6) 

" i 

It is customary to use the parameter, q — 1 — a. The above power-law tailed form of 
energy distribution can be fitted to experimentally observed particle spectra, and this way 



a numerical value for the parameter a can be obtained. The a = (q = 1) case recovers the 
classical Boltzmann-Gibbs-Shannon (BGS) formula 25N28] 

S B g = J2-Wi\nwi. (7) 

i 

According to this the quantity L(S) is to be maximized with constraints. Identifying 
the analogous formal logarithm for leading order non-additive energy composition, E12 = 
E 1 + E 2 + aE x E 2 , as 



one considers 



L(E) = - ln(l + aE), (8) 
d 



— In ^2 w] a — (3 ^2 w i~ l n (l + a Ei) — a^^Wi = max. (9) 



The maximum is achieved by the canonical probability distribution 

Wi = A(b(a + PLi))-* (10) 

with 

L i = -hi(l + aE i ), A = e-' L{ - S \ b = -^. (11) 
a 1 — a 

Then the normalization, the average and the definition of the entropy lead to the condition 

1 = ba + b(3 (L) . (12) 
Finally the equilibrium distribution simplifies to 

Wi = \(l + a(3Ly lla (13) 



Z 

with Lj given in eq. tfTTj) . Here we have introduced the following shorthand notations: 

Z = K-a(W» - <14) 

We should keep in mind that the reciprocal temperature, distinguished by the Zeroth Law, 
is the Lagrange multiplier (3. This is reflected well by the whole formalism, because the 
usual thermodynamic relations are valid. 

It is particularly interesting to consider now cases, when only one of the two quantities 
is composed by non-additive rules. In the limit of additive entropy but non-additive energy 
(a — > 0) the canonical distribution approaches 

^ = ^-(1 + 0^)-^°, where In Z = S BG - f3 (E) . (15) 
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Here Sbg is the Boltzmann-Gibbs-Shannon entropy (cf. eqJT]). For non-additive entropy and 
additive energy on the other hand a similar, but differently parametrized power-law tailed 
distribution emerges: 

^ = ^(l + /3a£ i )" lA \ (16) 

with 

= l-a(l + /?(£»• (17) 
The latter relation can be transformed into a more suggestive form by using q = 1 — a and 
the temperature parameters T = 1/(5 and T = 1//3: 

T=lf+0-l)(E). (18) 

By using the distribution given in eq. flIB"]) . the expectation value of the energy, (E), is directly 
given as a function of T and a = 1 — q. 

III. SUPERSTATISTICAL MONTE CARLO METHOD 

In either case discussed in the previous section, the generalized canonical distribution of 
the different energy states in a system in thermal equilibrium with non-additive composition 
rules is given by a formula 

m , = w 1+ ^r (19) 

Z T s \ c J 

In the c — > oo limit this formula coincides with the familiar Gibbs factor: 

lim Wl = -±-exp(-PEi). (20) 

The quantity q — 1 — 1/cis called the Tsallis index. Here c = /3/a and is in fact the inverse 
absolute temperature for the energy non-additivity case; for the entropy non-additivity on 
the other hand (3 has to be replaced by (3 and c by 1/a as it was explained in the previous 
section. The thermodynamic temperature in the latter case, according to the Zeroth Law, 
can be obtained by using eq.flTBl. 

The Tsallis distribution weight factor, Wi, on the other hand can be obtained as an 



integral of Gibbs factors over the Gamma distribution [29|, |30] , 

1 r°° 

Wi = — dOw c (e)exv{-epEi), (21) 

ZiTSt Jo 
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with 



w c (B) = ^—9^ e" c6 '. (22) 



T(c) = (c — 1)! for integer c is Euler's Gamma function. By its definition the integral of 



is normalized to one. This approach is a particular case of the so called super statistics 



Based on this, any canonical Gibbs expectation value, if known as a function of /3, can be 
converted into the corresponding expectation values with the power-law tailed canonical en- 
ergy distribution. The respective partition functions, Zq and Z?s ensure the normalization 
of the Wi probabilities, Y^,i w i — 1- They are related to each other: 

poo poo 

ZtsW) = J2 d6w c (6)exp(-6f3E i )= d9 w c (6)Z G (6f3). (23) 
i Jo Jo 

The above formula can be interpreted as averaging over different #/3-valued Gibbs simu- 
lations. The averaging is understood in the partition sum, meaning that the weighting 
'Boltzmann'-factor is also fluctuating. It assumes that the underlying process of mixing 
different inverse temperatures is much faster than the averaging itself. 

The question arises, which strategy is the best to follow in order to perform lattice field 
theory simulations with power-law tailed statistics instead of the Gibbs one. Neither the 
ensemble of different f3 values (Euclidean timelike lattice sizes), nor the re-sampling of the 
traditional, Gibbs distributed configurations is practicable in a naive way. The N t lattice 
sizes are limited to a small number of integer values - hence the good coverage of a Gamma 
distribution with an arbitrary real c value is questionable. The already produced configura- 
tion ensembles were selected by a Monte Carlo process according to the Gibbs distribution 
with the original lattice action; there is no guarantee that the re-weighting procedure (which 
includes part of the weight factors in the operator expressions for observables) is really con- 
vergent (i.e. does not contain parts growing exponentially or worse). We choose another 
strategy: we use 9 values selceted as random deviates from an Euler-Gamma distribution 
during the Monte Carlo statistics. 

The lattice simulation incorporates the physical temperature by the period length in the 
Euclidean time direction: (3 = N t at- Due to the restriction to a few integer values of Nt, we 
simulate the Gamma distribution of the physical — 1/T values by a Gamma distribution 
of the timelike link lengths, a t . We assume that its mean value is equal to the spacelike 
lattice spacing, a s . Then the ratio 9 = a t /a s follows a normalized Gamma distribution with 



the mean value 1 and a width of 1/ 1 \fc. (In the view of ZEUS e + e~ data c ~ 5.8 ± 0.5, the 
width is about 40 per cent.) In our numerical calculations we apply the value c = 5.5. 

For calculating expectation values in field theory a generating functional based on the 
Legendre transform of Z is used. Our starting assumption is the formula ( 1231) with 

Z G [9f3] = JvUe- s[u ' e] . (24) 

Since we simulate the canonical power-law distribution by a lattice with fluctuating asym- 
metry ratio, there are two limiting strategies to execute the Legendre transformation: i) in 
the annealing scenario the lattice fluctuates slowly and one considers first summations over 
field configurations, in the ii) quenched scenario on the contrary, the lattice fluctuations are 
fast, form an effective action (virtually re- weighting the occurrence probability of a field 
configuration), and the summation over possible field configuration is the slower process 
performing the second (i.e. the path-) integral. In this paper we investigate numerically 
the general case when one may choose when a new value for 9 is taken. The frequency of 
these fluctuations may go from one in each Metropolis step for the field configurations to 
one in the whole Monte Carlo process (the latter being the traditional method). Our results 
presented in the next section belong to a choice of 5 field updates for the whole lattice before 
choosing a new 9. This peculiar value was controlled by a series of simulations and proved 



to be sufficient for a close equilibration to a given, momentary temperature [31 ]. 

The effect of 9 fluctuation is an effective weight for field configurations, which may depend 
on a scaling power according to the time (or energy) dimension of the operator under study. 
In general we consider the Tsallis expectation value of an observable A[U] over lattice field 
configurations U. A may include the timelike link length, say with the power v. A — 9 V A. 
The Tsallis expectation value then is an average over all possible at link lengths according 
to a Gamma distribution of 9 = a t /a s . We obtain: 

(A)ts = — 7—7- [ d9 9 c ~ l e- c 9 [ VUA [U] 9 v e - s[e ' u] (25) 

z TS r(c) J J 

with 

Zts 



Jd9 9 c - 1 e- ce JvUe-WM. (26) 



T(c) 

The 9 dependence of the lattice gauge action is known for long: due to the time derivatives of 
vector potential in the expression of electric fields, the "kinetic" part scales like a t al/(afa1) = 
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a s /a t , and the magnetic ("potential") part like a t a? s / \a 2 s a 2 s ) = a t /a s [34]. This leads to the 
following expression for the general lattice action: 

S [9, U\ = a 9 + b/9, (27) 

here a = S SS [U] contains space-space oriented plaquettes and b = S ts [U] contains time- 
space oriented plaquettes. The simulation runs in lattice units anyway, so actually the U 
configurations are selected according to weights containing a and b. In the c — > oo limit the 
scaled Gamma distribution approximates 5(9 — 1), (its width narrows extremely, while its 
integral is normalized to one), and one gets back the traditional lattice action S = a + b, 
and the traditional averages. For finite c, one can exchange the 9 integration and the 
configuration sum (path integral) and obtains exactly the power-law-weighted expression. 

IV. STATISTICS OF POLYAKOV LOOPS 

Before discussing our results for the SU(2) pure gauge lattice field simulation using Euler- 
Gamma distributed timelike lattice spacing (and simulating this way a fluctuating inverse 
temperature to leading order in non-extensive thermodynamics), let us present a figure 
about the numerical quality of this randomization. In Fig{TJthe evolution process and the 
frequency distribution of the 9 values are shown for the reference run with c = 1024.0 and for 
the investigated case with c = 5.5. We have choosen a new value for the asymmetry ratio 9 
in each 5-th Monte Carlo update - in order to leave some time for the relaxation of the field 
to its thermal state at each instantaneous (59 inverse temperature. In the figure only each 
5-th value is shown. The Monte Carlo simulations were done at the coupling 4/g 2 = 2.40 
for this particular statistics with the Metropolis method. 

Our reference case, thought to be close to the c = oo traditional system, is specified by 
c = 1024. The re-fit to the distribution of effectively used values after 20000 draws from the 
Euler-Gamma distribution by a numerical subroutine was done by the statistics tool "gretl" . 
In the special case of our random weighting one expects an Euler-Gamma distribution with 
reciprocial a = c and — 1/c parameters. On the basis of a sample of 20000 9 values we 
achieved a reconstruction of a = 1009.8 and 1//3 = 1010.1. Similarly for c = 5.5 we obtained 
a = 5.5179 and 1//3 = 5.5255. 

Now let us turn to the discussion of the behavior of the order parameter of the confinement 
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FIG. 1: The Monte Carlo evolution and distribution of t asym = a t /a s for the coupling 4/g 2 = 2.40. 
The random deviates for the results shown in the upper row are thrown with the parameters 
a = c = 1024.0 and = 1/c = 0.000977. The re-fit by gretl gave a = 1009.8 and /3 = 0.000990. 
The same parameters in the lower row are a = c = 5.5 and /3 = 1/c = 0.181818. The re-fit by gretl 
gave a = 5.5179 and (3 = 0.18098. 

- deconfinement phase transition. The Polyakov Loop is calculated by taking the trace of 
the product of gauge group elements on timelike links closing a loop due to the periodic 
boundary condition: 

N t 

P(x) = Trl[U t (t,x). (28) 
t=i 

The traditional order parameter of the phase transition is the expectation value of the volume 
averages for each lattice field-configuration during the Monte Carlo process. For the gauge 
group £77(2) this quantity is real: 

Ke? = SRe^?(i). (29) 

X 

In our present investigations the characteristic width parameter of 1/T-fluctuations is c = 
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5.5, corresponding to a relative width of 1/y/c ~ 0.43. As a reference the c = 1024.0 case is 
taken - here the relative width is about 1/a/c = 1/32 w 0.03. 

The plots in Figj2] show the fluctuations of the order parameter 3fte P for the reference 
runs with c = 1024.0. The fluctuating values as a function of the Monte Carlo step are 
plotted on the left hand side, while their probability distributions on the right hand side. 
The values for the inverse coupling include both the confinement and deconfinement phases. 

By producing these results we took five consecutive Metropolis sweeps over the whole 
4-dimensional 10 3 x 2 lattice while keeping the asymmetry value 9 = a t /a s constant. Then a 
new 9 was chosen as a random deviate from an Euler-Gamma distribution. Only these 5-th 
values are plotted and counted for obtaining expectation values. The probability distribu- 
tions of these values were determined by using the statistics software tool "gretl'. Hereby the 
first 5000 configurations were sometimes taken out from the samples, consisting of 100000 
lattice configurations each, this did not change expectation values appreciably. For the 
statistical evaluation only each 5-th configuration was selected, being fairly independent of 
each other in the evolution governed by the Metropolis algorithm and certainly belonging to 
different 9 values. The frequency distributions reflect cleanly when several 9fte P expectation 
values are occurring during the Monte Carlo evolution, by several maxima. This is the case 
near to the phase transition point. 

Similar pictures from Monte Carlo simulations with fluctuating inverse temperature using 
the parameter c = 5.5 are plotted in the figures |3]-[7l Here the effect of the width in the 
possible temperature values is clearly seen in the larger fluctuations of the order parameter 
compared to the reference case c = 1024.0 at the same coupling. Also the critical inverse 
coupling strength moves towards higher values for c = 5.5. In Figj5]we zoom to the neigh- 
borhood of the critical coupling: The distribution of the 3fte P values are characteristically 
wide. In the third row, at 4/ ' g 2 = 2.14, the distribution of possible values is almost flat be- 
tween —1 and 1. (Due to the SU(2) trace normalization, as we use it, the maximal absolute 
value of the order parameter is 2.) The intermittent behavior between positive and negative 
values of JieP, a sure sign of the restoration of the center symmetry Z2, can be catched 
until the value 4/g 2 = 2.20, as it can be inspected in Figj6) For even higher inverse coupling 
strength the observational sample is too short to observe this effect. 

How to estimate the critical coupling for the appearence of the nonzero order parameter? 
The method closest to the traditional one 32] is to take the average value over the statis- 
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FIG. 2: The Monte Carlo evolution and distribution of IfceP for the couplings 4/g 2 = 
1.80, 1.85, 1.90 and 1.95 using c = 1024.0 from the top to the bottom. This reference pictures 
show a nearly-traditional confinement - deconfinement phase transition for the SU(2) Yang-Mills 
system. Note the small width of the order parameter distribution. 
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FIG. 3: The Monte Carlo evolution and distribution of IfceP for the couplings 4/g 2 = 
1.80, 1.85, 1.90 and 1.95 using c = 5.5 from the top to the bottom. Confinement phase. 
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FIG. 4: The Monte Carlo evolution and distribution of IfceP for the couplings 4/g 2 = 
2.00, 2.05, 2.06 and 2.08i using c = 5.5 from the top to the bottom. These couplings are nearly 
critical. 
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FIG. 5: The Monte Carlo evolution and distribution of IfceP for the couplings 4/g 2 = 
2.10, 2.12, 2.14 and 2.15 using c = 5.5 from the top to the bottom. Here the two-peak distribution 
develops, the deconfinement sets in. 
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FIG. 6: The Monte Carlo evolution and distribution of IfceP for the couplings 4/g 2 = 
2.16,2.18,2.20 and 2.25 using c = 5.5 from the top to the bottom. By these couplings we dwell 
into the deconfinement regime. 
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FIG. 7: The Monte Carlo evolution and distribution of IfceP for the couplings 4/g 2 
2.40, 2.45, 2.50 and 2.55 using c = 5.5 from the top to the bottom. For these couplings o 
one symmetry breaking maximum occurs representing a well-developed deconfinement phase. 
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FIG. 8: Results on Polyakov Loop spatial average expectation values in long runs (100.000 Monte 
Carlo steps, each 5-th kept) on 10 3 x 2 lattices at c = 5.5 (red squares) and at c = 1024.0 (green 
circles). The Gaussian widths are indicated by error bars. The transition point, i.e. the critical 
coupling strenth, x = 4/g 2 , is estimated by a functional fit, JJeP ~ (x — x c ) 1//3 . 

tics. In Figj8]we plot (dteP) over the longer Monte Carlo runs presented above with their 
distribution. There is a characteristic difference between the c = 5.5 and the c = 1024.0 
cases. A possible fit to the average values is given by a fractional power; it seems that a 1/3 
power-law behavior describes the critical scaling well. Of course, on the basis of the present 
data a square root behaviour also cannot be excluded. The obtained positions of the critical 
couplings differ: A/gl « 1.85 for c = 1024.0 while 4/g 2 m 2.12 for c = 5.5. 

For drawing conlcusions relevant to the physics the inverse lattice couplings have to be 
related to temperatures. Figure [TO] presents T/T c ratios versus the inverse coupling, 4/ q 2 
for N t = 2 lattices, based on data for critical couplings on different iVj-sized lattices 33 ]. 
Although those simulations were carried out without temperature fluctuations, i.e. taking 
c = oo, we use them as a first estimate for the temperature - coupling correspondence. The 
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4-th order cumulant for Polyakov Loop (10 3 x 2) 




FIG. 9: Fourth order cumulants of the Polyakov Loop spatial average in long runs (100.000 Monte 
Carlo steps, each 5-th kept) on 10 3 x 2 lattices at c = 5.5 (full circles) and at c = 1024.0 (open 
circles). The critical coupling strenth, x = 4/^, is obtained by a linear fit to the smaller nonzero 
values. 

critical coupling in our calculation for c = 1024.0 is close to the result obtained previously 
on same sized (N t = 2) lattices. The critical coupling at c = 5.5 - following the c = oo line 
of constant physics - corresponds on the other hand to a temperature which is 1.3 times 
higher than the usual value. 

V. CONCLUSION 

1. For c = 5.5 (a realistic value from p T spectra) the critical coupling at the deconfme- 
ment phase transition shifts towards higher values. To this value an increase of the 
deconfmement temperature is obtained at T c (5.5) 1.3T C (1024) 1.3T c (c = oo). 

2. Aiming at the same 1/T value for the simulation, i.e. (6) = 1, the temperature is 
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SU2 temperature vs inverse coupling 
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FIG. 10: Deconfinement temperatures based on [33J] vs inverse coupling strength obtained from 
different size lattice simulations (Nt values are indicated on the plot). The arrows point to our 
findings of critical temperatures with Nt = 2 for c = 1024.0 and c = 5.5, respectively. The 
corresponding horizontal lines are drawn at 1.00 and 1.30 with respect to the c = oo case. 

expected to make an increase of about 20 per cent due to (1/0) = c/(c — 1) ~ 1.22. 
This shows the same trend as obtained by the Monte Carlo simulations, but not its 
whole magnitude. 

3. We obtained, assuming the traditional scaling dependence between coupling and phys- 
ical temperature, an increase of 15 per cent in leading to about an increase of 
30 per cent in T c . The dynamical effect is definitely larger than the trivial statistical 
factor of 1.22. 

4. Therefore experiments aiming at producing quark matter under circumstances char- 
acteristric to high energy collisions should consider the possibility of an about 30 per 
cent higher T c then predicted by traditional Monte Carlo lattice calculations. A pos- 
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sible measurement of the value of the width parameter c can be achieved by analyzing 
event-by-event spectra. 

These preliminary conclusions are based on a comparison with the c = oo traditional 
results. In future works we aim to explore the T/T c — 4/g 2 curve and possibly the renor- 
malization of physical quantities under the condition of fluctuating temperature with finite 
c values. 
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